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We propose a new variational Monte Carlo (VMC) method with an energy variance extrapolation 
for large-scale shell-model calculations. This variational Monte Carlo is a stochastic optimization 
method with a projected correlated condensed pair state as a trial wave function, and is formulated 
with the M-scheme representation of projection operators, the Pfaffian and the Markov-chain Monte 
Carlo (MCMC). Using this method, we can stochastically calculate approximated yrast energies and 
electro-magnetic transition strengths. Furthermore, by combining this VMC method with energy 
variance extrapolation, we can estimate exact shell-model energies. 

PACS numbers: 21.60.Cs, 21.60.Ka 



Shell-model calculations have been a central issue in 
studies of nuclear structure, and constant effort has been 
made to solve large-scale shell-model problems. Exact di- 
agonalization, which is a standard method, has recently 
been able to handle large-scale problems with O(10 10 ) 
dimension [J However, computational feasibility in 
this method is limited and strongly depends on the size 
of single particle space and valence nucleon numbers. To 
overcome this problem and to extend the feasibility of 
shell-model calculations, various methods [§-[il| with dif- 
ferent kinds of algorithms have been proposed and im- 
proved. 

In this paper, we propose a new method for shell- 
model calculations using the Markov-chain Monte Carlo 
(MCMC). This new method is a variational Monte Carlo 
(VMC) with a projected correlated condensed pair state 
as a trial wave function. This new VMC can stochas- 
tically give not only energies but also electro-magnetic 
transition strengths, although these values are approxi- 
mated. To estimate the exact shell-model energies over 
the limitations of variational formulation, we use the en- 
ergy variance extrapolation, which has been studied in 
Refs.[H, [HI and has been successfully applied to the nu- 
clear shell model 0, [TTJ] . As this extrapolation needs a 
series of systematically approximated wave functions, we 
introduce a truncation scheme based on the spherical ba- 
sis in the form of a projection operator. 

In this study, projection operators are implemented by 
using the Monte Carlo in a novel way. In fact, particle 
number projection, magnetic quantum number projec- 
tion, parity projection, and projection onto truncation 
spaces can be implicitly performed without numerical 
integrations, and angular momentum projection is per- 
formed by two-dimensional numerical integration. 

Here we consider a new variational formulation of shell- 
model calculations. As a trial wave function for a nuclei 
with N v valence protons and N v valence neutrons, i.e., 



N = + N u , we take as 

= GP\4>) 

and 
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where / is a skew-symmetric matrix. The |0) is an inert 
core and the cj 's are proton or neutron creation operators 
{i = 1, • • • , Nir : proton and i = N n + 1, • • • , N : neu- 
tron). By this parametrization, proton- neutron pairing 
correlation is included, in addition to the proton-proton 
and neutron-neutron pairing correlations. For simplicity, 
to present this formulation, we assume that the compo- 
nents of / are real numbers. We can also use a complex 
number although this would require some modifications. 
The P is a projection operator and we use the following 
two kinds. One is 



P = P Io P*P, 



M 



(3) 



where P Ia , P v and Pm are projectors of the ^-component 
of isospin Iq, parity tt and ^-component of angular mo- 
mentum M, respectively. The other is 



p _ plo pi* pJ 
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where P J is a projection onto angular momentum J. The 
G is a correlation factor as 



(5) 



where cv's are variational parameters. The rii is taken as a 
number operator for orbit i, so that the G is commutable 
with the angular-momentum projector. The G can be 
easily evaluated as will be discussed later. 

For the trial wave function in Eq.Q, the energy ex- 
pectation value is given by 



E=(H) 
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To evaluate this, we introduce an alternative representa- 
tion of the projection operator in Eq.Q, which can be 
shown by the so-called M-scheme states defined as 



to) 



.10), 
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where to = (mi, W,2, ■ • ■ m^). As the M-scheme state 
is an eigenstate of the z-components of the isospin and 
angular momentum, and parity, the projection operator 
is given by 



P Io P*P M = 



E 
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where the summation is restricted within the M-scheme 
states with given quantum numbers, that is, particle 
numbers N„ and N u , and magnetic and parity quantum 
numbers M and tt. By this representation, {m\P\<j)) can 
easily be calculated as 



{m\P Ia P^P M \<ty) = (m\4>)5 z(m)>M , 



(9) 



where z(m) gives the parity and magnetic quantum num- 
bers of the M-scheme state |m). 

By introducing the M-scheme representation of the 
projection operator, the G-factor becomes a c-number 
because the G-factor is diagonal in the M-scheme, that 
is, 



G|m) = G{m)\m), 



(10) 



where the G(m) is an expectation value of the G con- 
cerning the M-scheme state m). Moreover the matrix 
elements of the H between M-scheme states can also be 
evaluated, and we denote them as 



(m\H to') = h ri 



(11) 



where to) and \m') have the same quantum numbers and 
h is generally very sparse because a shell-model Hamil- 
tonian consists of one-body and two-body interactions. 

By these relations (| 10[) and (fTTj) , Eq.© can be rewrit- 
ten as 



E 



p(m)E L (m), 



where the local energy _E/,(m) is defined as 

(m'\P\<f>)G(m') 
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and the sampling density p(m) is defined by 
\(m\P\4>)G(m)\ 2 
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where p(m) > and Em£{M"} P( m ) = l - 

In numerical calculations, the energy formula Eq. (|12[) 
is not practical because the dimension of the M-scheme 



space becomes intractably huge as the size of single par- 
ticle space and proton and neutron numbers increase. 
Eq. ([T2|) . however, is suited to Monte Carlo calculations 
because if we can generate a set of the M-scheme states 
|to) that obey the occurrence ratio p(m) by the Monte 
Carlo sampling, the energy expectation value can be es- 
timated by 



(15) 



where the No is a number of Monte Carlo samples. 

In this formulation, the projection operator appears 
only in the projected overlap between the and the M- 
scheme state \m). In the case of the projection operator 
Eq.([3]), the projected overlap simply becomes the overlap 
(to \4>) as shown in Eq.([9]), and this overlap can be given 
by the Pfaffian as 



N 



(16) 



where (X 
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[15j. The definition of the 
Pfaffian is shown in the Appendix. 

Next we delve into the Markov-chain Monte Carlo. To 
evaluate Eq. (IT5|) by the Monte Carlo method, a set of \m) 
whose distribution obeys Eq.([TJ} is needed. Here we con- 
sider a random walker \m) on the M-scheme space with 
the given proton and neutron numbers, parity, and mag- 
netic quantum numbers. A random walker |to) moves 
to \m!) on the M-scheme space with the same quantum 
numbers in the following way: 

• We choose two nucleons (m,,i7ij) in the \m) ran- 
domly and annihilate the nucleons in the |to). We 
call the resultant state |toj) = c mj c mi \m). 

• We sum up the magnetic quantum numbers, the 
z-component of the isospin and the parity of the 
chosen nucleons. 

• In | to,/), we randomly choose two available unoc- 
cupied states, (m'^m'j), whose summed quantum 
numbers are the same as those of (to^to^). Then, 
we create two nucleons on chosen states. We call 
the resultant state Ito') 



,<r , \mj) 



Transition of random walker can be controlled in two 
ways. One is the Metropolis-Hasting (MH) algorithm. 
Whether or not a random walker \m) moves to \m'} de- 
pends on the ratio p(m / ) as 



p(to') 



(m'\P\^)G(m') 
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If p(m') > 1, the walker \m) always moves to \m'). If 
p(m') < 1, according to the p(m'), we determine whether 
or not the walker \m) moves to \m'}. 

The other is the Gibbs sampling algorithm. When we 
consider a random walker |m), we choose a removed pair 
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randomly and obtain |mj). Then, we calculate allp(ra')'s 
for possible m') = cr ,cl ,\mj) conserving appropriate 

i j 

quantum numbers. Finally, we choose \m') by transition 
probability W(mi — s- m') defined by 



W(mi — ¥ to') 



p(m') 
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These algorithms satisfy detailed balance and ergodicity. 

Next we consider how to optimize the parameters' a's 
and /'s in the trial wave function. To use the steepest 
descent method, a gradient vector needs to be evaluated 
by the Monte Carlo method. Each component of the 
gradient vector is given by straightforward calculations 
as 



dE 

dan 
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where 0/ 4 . is an operator, the matrix elements of which 
are as follows 



O hj {m) 
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Sophisticated derivation is shown in Ref . [la ] . The deriva- 
tion needs some extensions if we use complex numbers for 
fij. The gradient vector obtained by the Monte Carlo 
method suffers from stochastic noises. To reduce such 
noises we use the stochastic reconfiguration (SR) method 
14 1, the details of which are shown also in Ref. [l5| . In 
this way, we can numerically evaluate the gradient vec- 
tor and can optimize the parameters of the present wave 
function based on the steepest descent method. 

In this formulation, proton and neutron number pro- 
jections, magnetic quantum number projection, and par- 
ity projection can be implicitly performed without nu- 
merical integration, while angular momentum projection 
can not. The J-projection Pj^ [Tt| is given by 



P J 
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where the f2 stands for Euler's angles (a, (3, 7) and J — 
2J+1. The i?(J7) is a rotational operator and D J MK {Vt) is 
Wigner's /^-function. The g's are additional variational 
parameters. 



The angular momentum operator is rewritten by 



(25) 



where 



pJ 



4tt 
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The d J MK (fi) is Wigner's d- function. The projected over- 
lap {m\P l0 P* P(j\4>) becomes 

{m\P I °P'Pi t \4>) ^Y.9K{MP J MK\<t>)5 z[m)M ^ (27) 

K 

with Pm = YIk 9 k ^M k- The projected overlap con- 
cerning P(j K is 



{m\P J MK \ 
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with 



d$dl sin pd J MK (f3)e- lK ^R(/3, 7) 

(28) 



R(j3,j) = (m\e iJ y^e iJ ^\ 



(29) 



which can be calculated by the Pfaffian. Thus, the angu- 
lar momentum projection can be performed by the two- 
dimensional integration, which is a distinguishable fea- 
ture in this formulation. 

Finally, we consider how to compute electro-magnetic 
transition strengths. To evaluate them, angular momen- 
tum projection is indispensable and unnormalized initial 
and final states can be denoted by 

|W; J*,M a ) = GY,9K.Pfc,, K J<f>,) (30) 

where a = i, f. The Ji and Jf are spins of the initial and 
final states and the parameters of these wave functions 
are f's, a's and g's. 

The B(E2; Ji — > J/), for instance, is defined as 



B(E2; Ji ->• J f 



1 



2J t . + 1 
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where Q is a quadrupole operator, and \ipi) and \ipf) are 
normalized initial and final wave functions with spins Ji 
and Jf, respectively. For this Monte Carlo evaluation, 
noting the following relation as 



$f\Q\4>i) 



(32) 



where \ij) a ) = \ipa) / \f (i^alipa), we can perform the 
and ^jj^rj^y-, respectively. For the 



MCMC for 



MCMC of the {j $j!jry term with Eq.(|30| . the sampling 
density p(jn) is given by 



p{m) 
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and local quadrupole strength Q_L(m) is given by 
Qi(m) = 



nt- , {m'\Ph \6i)G(m') 



m'£{M 4 '} 



where =Ejr B ^Pj 



(m|P^|0/)G(m) 



Note that we take differ- 



ent M-scheme spaces for |m) (J z |m) = M/|m)) and |m') 
(J z |m') = Mi\m')). Thus, the strength of the electro- 
magnetic transition can be evaluated in the MCMC. 

Next, we numerically investigate the present varia- 
tional Monte Carlo method. Here we consider the 56 Ni in 
the pf shell with the GXPF1A inter action pj]], of which 
dimensions are about 1.09 billion in the M-schcmc. 

First, we carry out the VMC with a J z = J space where 
we consider the state with angular momentum J. To per- 
form the VMC, we prepare an initial wave function ran- 
domly and simulate the sampling density p(m) in Eq.(fT4|) 
by the Monte Carlo. Here we use the Gibbs sampling al- 
gorithm. After appropriate burn-in steps (~ 1000), a 
random walker moves more than 5000 steps in the M- 
scheme space. These numbers depend on the acceptance 
ratio and required accuracy of numerical calculations. In 
this way, we prepare several tens random walkers and es- 
timate the energy and energy gradient vector with statis- 
tical errors. With the aid of the SR technique, we modify 
all the variational parameters of the wave function and 
repeat this optimization process until the energy varia- 
tion goes to zero. In Fig. [TJ we show the convergence 
patterns of the energies with J z = 0, 2 and 4 states as 
functions of the iteration number for 56 Ni. The statisti- 
cal error during the optimization procedure is a few tens 
keV, which is too small to be shown in this figure. 

Next, we stochastically calculate the J-projected en- 
ergy by carrying out a J-projection on the converged 
wave functions at J z = J space. We call this stochas- 
tic VBP (variation-before-projection). Note that the pa- 
rameters of the wave function are optimized concerning 
number, parity, and J z projected energy. Here, we take 
9K = Sk.j for simplicity because we carry out the J- 
projection onto the wave functions optimized in the space 
with J z = J. For J = state, after 1000 steps as a 
burn-in, a random walker moves 500000 steps in the M- 
scheme space. We evaluate the energy by 10 random 
walkers. The energy is -205.333 ± 0.004 MeV for 56 Ni. 
We show the J-projected energies for J = 0, 2, and 4 
with the label VBP in Fig. [TJ The J-projection improves 
while there are still sizable differences between these VBP 
energies and exact shell model energies. In this formu- 
lation, we can stochastically evaluate electro-magnetic 
transition strengths with statistical errors. The calcu- 
lated B(E2)'s with the VBP wave functions for -> 2 
and 2 — 4 are 690 ± 8 and 134 ± 2 e 2 /m 4 , respectively. 
Here we use effective charges e„ = 1.5 and e u = 0.5. 

To overcome the variational limitation of the VMC, 
we introduce the energy variance extrapolation, which is 
a technique to estimate the exact energy from a series 
of approximated wave functions in a well-controlled way. 
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FIG. 1: (Color online) Convergence patterns of energies with 
J z — 0, 2 and 4 states as functions of the iteration number for 
56 Ni in the pf shell with the GXPF1A [r| interaction. The J- 
projected energies (VBP) onto the converged states at J z = 
J space are shown with the extrapolated (Ext.) and exact 
(Exact) shell model energies. In the inset, extrapolations of 
energies with J = 0, 2 and 4 as functions of energy variance 
are shown. 



This technique is based on a well-defined scaling prop- 
erty for energy eigenvalues. We define a difference 5E 
between energy eigenvalue (H) in a given subspace and 
exact energy eigenvalue (H)q, that is, 6E = (H) — (H)q, 
and an energy variance AE in the subspace is also de- 
fined as AE = (H 2 ) — (H) 2 . The difference 6E vanishes 
linearly or quadratically as a function of the energy vari- 
ance AE [9J]. By this scaling property, we can estimate 
the exact shell-model energies as the limit of zero energy 
variance. 

To apply this technique to the VMC, we introduce a 
truncation scheme in the form of a projection operator. 
In the case of the pf shell, due to a relatively large gap 
of spherical single particle energy between the / 2 orbit 
and others r {f^/%, P3/2 an d P1/2 orbits), particle- hole 
excitations across this shell gap form truncation spaces, 
®s<t(f7/2) A ~ 40 ~ s (r) s - The P* is a projection onto this 
truncation space and is added to the projection in Eqs.([21 
H]). This projection operator is quite simply realized 
by the restriction of the summation in Eq.® and we 
can easily include this truncation scheme in the present 
VMC. 

It is noteworthy that, in the Monte Carlo calculations, 
we find a fast computational method of the energy vari- 
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ance within the truncated space t as 
ie(t+2) 

where £x is the local energy defined in Eq. flT3l) . ATq is a 
number of Monte Carlo sampler, and Wt is a reweighing 
factor defined by the ratio between the number of random 
walkers within the (t + 2)-space, that is, No, and one of 
random walkers within the i-space. By this sampling 
method, we can avoid the explicit calculation of four- 
body interaction for energy variance and its computation 
becomes possible. 

In the inset of Fig. [TJ energies of the VMC calculations 
with different truncation spaces and different initial con- 
ditions are plotted as functions of the energy variance, 
AE. To the limit of zero energy variance, we can linearly 
extrapolate the energies with statistical errors obtained 
by the \ 2 fitting. The extrapolated energies agree with 
the exact ones. In the same way, we can extrapolate the 
B(E2). 

In conclusion, we have proposed a new variational 
Monte Carlo method with energy variance extrapola- 
tion for large-scale shell-model calculations. We have 
presented a formulation of wave function optimization 
based on the MCMC. In this method we can calculate 
approximated energy, other matrix elements, and electro- 
magnetic transitions for yrast states. Combining the 
VMC with energy variance extrapolation, we can esti- 
mate exact shell-model energies. This is an alternative 
extension of the VMC and is free from the sign-problem. 
By taking 56 Ni in the pf shell, we have shown the fea- 
sibility of large-scale shell-model calculations. Note that 
the present calculations can be carried out with a single 
core of the common PC. For larger computations, parallel 



computation plays a significant role. 

For further improvement of the present method, the 
stochastic VAP (variation-after-projection) concerning 
J-projection can improve accuracy of energy and energy 
variance to enhance the reliability of energy variance ex- 
trapolation. We are pursuing the stochastic VAP includ- 
ing its parallel computation as especially fitted to a state- 
of-the-art massive parallel computer, the results of which 
will be presented elsewhere. 

For energy variance extrapolation, we introduced the 
particle-hole truncation scheme into the VMC while we 
can use the arbitrary basis-truncation scheme. The 
excitation-energy truncation scheme in the No-core Shell 
Model [18] is also promising. We will investigate this 
direction in the future. 

Finally, this study is motivated by a study of the Hub- 
bard model[l5|. The present study sheds light onto a 
new way of projection in these stochastic calculations, 
an aspect of which may also be useful in applications of 
condensed matter physics. 

One of the authors (N.S.) was supported by Grants-in- 
Aid for Young Scientists (20740127) from JSPS and the 
HPCI Strategic Program from MEXT. 

Appendix 

The Pfaffian is defined for a 2n x 2n skew-symmetric 
matrix A as 

1 " 
P/(A) =2^ Z sgn^II^-iM^) ( 36 ) 

' ff£S 2 „ i=l 

where a is a permutation of {1, 2, 3, • • • , 2n}, sgn(cr) is its 
sign, Sm is symmetry group and et's are matrix elements 
of A. An efficient computation of the Pfaffian can be seen 
e.g. in Ref.[l|. 
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